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Abstract 

An overview of theories related to vibrational energy relaxation (VER) in proteins is pre- 
sented. VER of a selected mode in cytochrome c is studied using two theoretical approaches. 
One is the equilibrium simulation approach with quantum correction factors, and the other 
is the reduced model approach which describes the protein as an ensemble of normal modes 
interacting through nonlinear coupling elements. Both methods result in estimates of the 
VER time (sub ps) for a CD stretching mode in the protein at room temperature. The 
theoretical predictions are in accord with the experimental data of Romesberg's group. A 
perspective on future directions for the detailed study of time scales and mechanisms for 
VER in proteins is presented. 

Classification: Physical Science, Biophysics. 

Abbreviations: VER, vibrational energy relaxation; LTZ, Landau- TcUer-Zwanzig; Mb, 
myoglobin; cyt c, cytochrome c; QCF, quantum correction factor. 

1 Introduction 

When a protein is excited by ligand binding, ATP attachment, or laser excitation there occurs 
vibrational energy relaxation (VER). Energy initially "injected" into a localized region flows 
to the rest of the protein and surrounding solvent. VER in large molecules (including pro- 
teins) itself is an important problem for chemical physics |^ [2] . Even more significant is the 
challenge to relate VER to fundamental reaction processes, such as a conformational change or 
electron transfer of a protein, associated with protein function. The development of a accurate 
understanding of VER in proteins is an essential step toward the goal of controlling protein 
dynamics [Sj. 

Due to the advance of laser technology, there have been many experimental studies of VER 
in proteins @lEliniIZllHliniEl^Dll2,Ji3j J^^IEIEIEI- These experimental works are impressive 
but it is difficult to derive detailed information from the experimental data alone. Theoretical 
approaches including atomic-scale simulations can provide more detailed information. In turn, 
experimental data can be used to refine simulation methods and empirical force fields. This 
combination of experimental and theoretical studies of for protein structures and dynamics has 
begun to blossom. As experimental methods develop further and theoretical approaches grow 
in accuracy, the relationship will become fruitful. 

There have been many theoretical tools (Sec. [21) developed to analyze VER in proteins. 
Some aspects of VER in proteins can be explained by perturbative formulas based on the 
equilibrium condition of the bath (Sec. EJ, but the use of the perturbative formulas may be 
too restrictive to generally describe protein dynamics at room temperature. In this paper, we 

* fujisaki@bu.edu 
^ straub@bu.edu 



1 



future study of VER in proteins. 



2 Theories 

In this section, we present a selective overview of theories appropriate for the study of VER in 
proteins. For the most part, these theories have been developed to deal with VER in liquids, 
solids, or glasses. The reader is referred to a number of recent reviews |18| ll9 | On]. We refer to 
two distinct categories: one based on equilibrium dynamics and Fermi's golden rule, while the 
other is based on nonequilibrium dynamical models. 

2.1 Fermi's golden rule 

If (a) there is a clear separation between the system and bath, (b) the coupling between them 
is weak enough, and (c) the bath is assumed to be at thermal equilibrium, we can use quantum 
mechanical perturbation theory to derive a vibrational population relaxation rate through 
Fermi's golden rule [El EOj 



T{t) is the quantum mechanical force applied to the relaxing bond (system) considered, ms 
is the system mass, ujs is the system frequency, /? is an inverse temperature, and the above 
bracket indicates a quantum mechanical average. 

However, this time correlation function is very hard to numerically calculate. As a result, 
many approximate schemes have been proposed to address this limitation. A number of the 
most successful approaches is mentioned below. 

2.1.1 Landau- Teller-Zwanzig formula 

The most simple approximation is to take the classical limit {h 0) of Eq. 



Here the bracket denotes a classical ensemble average. This is called the Landau- Teller-Zwanzig 
(LTZ) formula, which has been applied to the study of VER in liquids [^. This strategy was 
used by Sagnella and Straub to discuss the VER of CO in Mb* CO [22] • This approxima- 
tion should be good for low frequency modes, but it becomes questionable for high frequency 
modes due to quantum effects. As such, advanced methods have been proposed to address this 
deficiency of the LTZ formula. 

2.1.2 Quantum correction factor 

The first alternative to the LTZ formula is the quantum correction factor (QCF) method. The 
basic idea of the QCF method is to relate a quantum mechanical correlation function with its 
classical analog |23| . When this is done for the force autocorrelation function in Eq. the 
final expression for the VER rate l/T^*^^ is 
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phonon relaxation process (harmonic QCF) 
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In the previous work ^lEOl, this result was expressed as T^'^^ ~ [(3huJs/Q{u!s)]Ti^ which is 
correct in the limit (3?iuJs ^ 1, as was appropriate for those studies. 

If the relaxation process is the linear resonance (1:1 Fermi resonance), then Q{ijJs) = 
Qh{^s)j i-G-) T^'^^ = Tf' [21]. Skinner and coworkers have provided a theoretical frame- 
work for organizing and expanding on a variety of QCFs appropriate for specific dynamical 
processes, dependent upon the underlying mechanism of VER. Though this strategy has been 
criticized |25 | 1261 l27j . it is known that the QCF method works rather well for specific problems 
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2.1.3 Reduced model approach 

An alternative approach to address the shortcomings of the LTZ formula is to use the re- 
duced model approach ^1^112111, which exploits a normal mode picture of the protein. By 
representing the Hamiltonian in terms of system, bath, and interaction terms 
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Calculating the force from this Hamiltonian, and substituting it into Fermi's golden rule Eq. 
we can derive a lowest order VER rate as |19[ 120] 
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and in previous papers HO], ms in the perturbative formulas should read ms = 1 as 
mass-weighted coordinates were employed. 
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to broaden the delta functions for numerical calculations. There exists another well known 
formula to describe the VER rate, the Maradudin-Fein formula |291 118j . 

W = TVdecay + VTcoll, (15) 
Wdecay = 7^ 2^ {l + nk+ni)^^— --^^ (16) 

VFcoU = — 2^ — {nk-ni)^^— ■ -2 (17) 

^ ^k^l 7 + (^^5 + ^k- ^lY 

with a width parameter 7. These two formulas are numerically similar for small 7, and equiv- 
alent for the limit of 7 — > 30^ This is a quantum mechanically exact treatment given the 
approximate truncated form of the interaction Hamiltonian. We have found that the trunca- 
tion error (the contribution from higher order terms) can be a serious problem, especially for 
proteins. For a more accurate treatment of VER, we must appeal to more advanced methods, 
described below. 



2.1.4 Other (advanced) approaches 

Methods that complement the above three methods involve calculating the force auto correla- 
tion function C,{t) appearing in Fermi's golden rule using different levels of approximations. Shi 
and Geva ,22, used a semiclassical approximation j^Jj for (^(t), and showed that even the slow 
relaxation of neat liquid oxygen (at 77K) can be well reproduced by their method. From their 
study, it was shown that the short time dynamics of C,{t) is important to predict the correct 
VER rate. This implies that the short time approximation may be adequate for an accurate 
estimate of C(t). Various time-dependent self-consistent field methods ^2] or path integral 
methods 33 should be applicable to calculate C,{t). For other methods, the reader is referred 
to additional works [Hil ESI IHEj • 

To derive Fermi's golden rule, we have used the Bader-Berne correction |^, which holds 
only for harmonic systems. Bader, Berne, Pollak, and Hanggi extended this to an anharmonic 
system within a classical framework |37j, and found that the VER of such a system can be 
nonexponential in time and is significantly affected by the character of the bath. This consid- 
eration will be important when one studies the VER of CO in Mb, especially for the VER of 
a highly excited CO bond. 



2.2 Nonequilibrium simulation 

The above equilibrium simulation methods based on Fermi's golden rule invoke several as- 
sumptions as described above. These assumptions might be invalid in some cases. As VER is 
a nonequilibrium phenomenon, the appeal of nonequilibrium approaches is quite natural. 



2.2.1 Classical approaches 

Classical nonequilibrium simulations to investigate VER in proteins were first conducted by 
Henry, Eaton and Hochstrasser 38 . In conjunction with their experimental studies, they 
employed classical molecular dynamics simulations of heme cooling in Mb and cyt c in vacuum 
and found that heme cooling occured on two time scales: short (1-4 ps) and long (20 ps for 
Mb and 40 ps for cyt c). Nagaoka and coworkers carried out the similar simulations for Mb in 
vacuum and obtained similar time scales 39 . Importantly, they found that the normal mode 
frequencies localized in the proprionate side chains of the heme are resonant with the water 
vibrational frequencies. 
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and Straub showed that the VER for Mb in water can be described by a single exponential 
with a few ps VER time |^ . Furthermore, they suggested that the main doorway of VER 
is due to the coupling between the proprionate side chains and water, which is in accord 
with Nagaoka's and Hochstrasser's observations. Bu and Straub supported this view through 
simulations of mutant Mb's and Mb variants having structurally modified heme groups 41_. 
They also investigated VER of cyt c in water, and found that the VER presents a biphasic 
exponential decay with two VER times: fast (a few ps) and slow (tens of ps) 

Kidera's group studied VER in proteins from a different perspective UHl- They excited a 
single normal mode in Mb, and examined the vibrational energy transfer (VET) between normal 
modes. As is well known, VET is caused by (nonlinear) Fermi resonance: if the frequency 
matching is good, and the coupling between normal modes is strong enough, there will be 
VET. This picture is very useful to characterize VET at low temperatures. However, at high 
temperature there occurs non-resonant VET. They numerically found that the amount of VET 
is proportional to a reduced model energy including up to third order coupling elements (see 
also SecEXHl). 

2.2.2 Quantum approaches 

For all but the simplest systems, quantum approaches for nonequilibrium simulations are ap- 
proximate and time-consuming. Nevertheless, these methods can overcome problems in inher- 
ent to classical simulations. There are two categories: vibrationally quantum methods, and 
electronically quantum ones. 

Hahn and Stock used a reduced model (consisting of the retinal rotation and other envi- 
ronmental degrees of freedom) to describe the pump-probe spectroscopy for the retinal chro- 
mophore in rhodopsin j44j . Flores and Batista, employing the same model, suggested the pos- 
sibility to control the retinal rotation by two (chirped) laser pulses '45'. To solve the quantum 
dynamics for the large system, they employed time-dependent self-consistent field (TD-SCF) 
methods |32j . Notably, vibrational SCF methods have been used to calculate vibrational energy 
levels for a small protein (BPTI) |46j . 

The combination of classical simulations for vibrational motions and quantum calculations 
for electronic structure, in some portion of the molecule, has been widely used for the cal- 
culations of up to moderate-sized molecules. One cutting edge application to a large system 
is the calculation of bacteriorhodopsin's photoisomerization in the excited chromophore state 
by Hayashi, Tajkhorshid, and Schulten 07j. In their treatment, a portion of the retinal chro- 
mophore including three double bonds was treated as the quantum mechanical region, and the 
complement, including the protein and water, as the molecular mechanical region. During the 
simulations, there occurs nonadiabatic transitions betweeen two electronic states (So and Si) 
which was treated semiclassically. They numerically showed that only one bond (Ci3=Ci4) 
rotates unidirectionally due to the coupling with the protein, and found that several other 
bonds can twist in any direction if there is no protein. 

3 Cytochrome c 

In this section, we shall focus on one protein, cytochrome c (cyt c), and review the recent 
theoretical studies about this protein. There are several reasons to select this protein as a 
prototypical one: (a) Cyt c is a relatively small protein with 1745 atoms. Other proteins 
of similar scale are Mb, BPTI, and human lysozyme. (b) The detailed X-ray structure is 
known for cyt c. (c) Cyt c has a function of electron transfer. The basic theoretical and 
computational works on cyt c were summarized by Wolynes and coworkers |48j . Wang, Wong, 
and Rabitz studied VER in cyt c using their hydro dynamical method j49j . Garcia and Hummer 
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describe the results of our studies on the VER of cyt c using two different methods (the QCF 
method in Sec. 12.1.21 and the reduced model approach in Sec. I2.1.3j) and compare them with 
the experimental results of Romesberg's group IT^ [T7| . 




Figure 1: Cytochrome c near heme. Only the 80th methionine (MetSO) residue and heme are 
depicted. Relevant atoms (C, D, S, Fe) are also indicated. This figure was created by VMD 
(Visual Molecular Dynamics) |51j . 

3.1 Quantum correction factor approach for cyt c 

Bu and Straub |S2] employed the QCF approach (Sec. I2.1.2j) to estimate the VER rate of a 
CD bond in the terminal methyl group of MetSO in cyt c (see Fig. Their calculations were 
carried out using the program CHARMM |53j . and cyt c was surrounded by water molecules 
at 300K. In Fig. |5J we show the force autocorrelation function and its power spectrum. With 
the CD bond frequency ujs = 2133 cm~^, we find = Qii^s) — 0.4 ~ 1.0 ps~^, so that the 

classical VER time is 1.0 ~ 2.5 ps. 
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Figure 2: Left: Averaged force auto correlation function for four trajectories at 300K. Right: 
Fourier spectrum for the four correlation functions with error bars. 

Since the CD bond frequency is located in a transparent region of the vibrational density of 
states, with no other state overlapping with this frequency [SI], it is concluded that there is no 
linear resonance (1:1 resonance). To use the QCF method, we thus need to assume nonlinear 
resonances corresponding to multiphonon VER processes. If the VER process assumes that 
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quantum of vibrational energy, the appropriate QCF (harmonic-harmonic QCF) is j23j 

QuHii-os) = QH{i^A)QH{'^s - wa)- (18) 

Alternatively, if the VER process is one that leads to the excitation of one bath vibrational mode 
of frequency uja^ with the remaining energy fi{u:s — ^a) being transferred to lower frequency 
bath rotational and translational modes, the appropriate QCF (harmonic-harmonic-Schofield 
QCF) is |2H] 

QH-Hs{i^s) = QHM^jQniiOs - a;A)e^^("^-"^)/^ (19) 

We need to determine the value of lva to use these formulas. From the normal mode and 
anharmonic coefficient calculations carried out in Sec. 13. 2( we have found that the CD mode is 
strongly resonant with two lower frequency modes, a;i655 = 685.48 cm~-^) and CJ3823 = 1443.54 
cm~^), where \iOs — wiess — ^^^3823 1 = 0.03 cm~^ for the standard parameters of CHARMM. 
We might be able to choose uja = 1443.54 cm^-*^ or 685.48 cm^-*^. If we choose uja = 1443.54 
cm~^ at 300K, Tf/T^'~'^ = Qhh{^s)/Qh{^s) = 2.3 for the harmonic-harmonic QCF and 
j.cyj.QCF _ Q^_^g(ijjgyQ^{ijjg"^ = 2.8 for the harmonic-harmonic-Schofield QCF. Thus we 
find T^^^ = TfV(2.3 ~ 2.8) ~ 0.3 ~ 1.0 ps. 

3.2 Reduced model approach for cyt c 

Fujisaki, Bu, and Straub ^| |2^ took the reduced model approach (Sec. I2.1.3|) to investigate 
the VER for the same CD bond stretching in cyt c. However, in their calculation, all modes 
represent normal modes, so the CD "bond" turned out to be the CD "mode." Using the 
formulas in Sec. 12.1.^ they calculated the VER rate for the CD mode {ucd = 2129.1 cm~^) 
and other low frequency modes (a'3330 = 1330.9 cm~^, wigge = 829.9 cm~^, wiess = 685.5 
cm~^) as a function of the width parameter 7 (Fig. inj. To this end, they needed to calculate 
anharmonic coupling coefficients according to the formula 

<^S,k,l - 7;^ — a — H~ - o UikUji — (2U) 

2 dqsdqkdqi 2 2Aqs 

where Uik is an orthogonal matrix that diagonalizes the (mass-weighted) hessian matrix at the 
mechanically stable structure Ky, and Kij{±Aqs) is a hessian matrix calculated at a shifted 
structure along the direction of a selected mode with a shift ±Aqs- 

If we take 7 ~ Auj ~ 3 cm~^, we have Ti ~ 0.1 ps, which agrees with the sub-picosecond 
time scale for relaxation predicted using the QCF method (T^^^^ = 0.3 ~ 1.0 ps). We also 
see that the low frequency modes have longer VER time, a few ps, which agrees with the 
similar calculations by Leitner's group |18j . In the right of Fig. we show the temperature 
dependence of the VER rate. At low temperatures, the VER rate becomes flat as a function 
of temperature. At these lower temperatures, the VER is caused by the remaining quantum 
fluctuation associated with zero point energy. 

3.3 Related experiment 

Here we discuss the related experiment by Romeberg's group ^1 El Ej • They measured the 
shifts and widths of the spectra for different forms of cyt c; the widths of the spectra (FWHM) 
were found to be AwpwHM — 6.0 13.0 cm~^. If we can neglect inhomogeneous effects, the 
estimate of the VER time becomes 

Ti ~ 5.3/A6JFWHM (ps) (21) 
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Figure 3: Left: VER rates for the CD mode {ujcd = 2129.1 cm^^) and the other lower 
frequency modes (W3330 = 1330.9 cm~^, ujiqqq = 829.9 cm~^, W1655 = 685.5 cm~^) as a function 
of 7 at 300K. Right: Temperature dependence of the VER rate for the four modes with 7 = 3 
cm~^. 



which corresponds to Ti ~ 0.4 ~ 0.9 ps. This estimate is similar to the QCF prediction using 
Eq. (0J) (~ 0.3 ~ 1.0 ps) and larger than the estimate by the reduced model approach using 
Eq. (|11|) or (|T5|) (~ 0.1 ps). This might be due to the strong resonance between the three 
modes (4357, 3823, 1655), which forms a peak near 7 ~ 0.03 cm~^ in the left of Fig. |31 This 
resonance causes an increase in the VER rate, so we can say that this estimate of the VER 
rate is too large. On the other hand, there is no peak for the low frequency modes for 7 < 10 
cm~^; the estimate of the VER rate does not seem to be affected by the resonances. 

Note also that Romesberg's group studied Met80-3D, methionine with three deuteriums on 
the terminal methyl group, while we have examined Met80-1D, with one deuterium. It is known 
that the CHARMM force field calculation does not give an accurate value of the absorption 
peak. On the other hand, the DFT calculation for the methionine leads to much better results 
(Matt Cremeens, private communication). Clearly, we must improve our force field parameters 
according to DFT calculations, and examine how further optimization of the parameters affects 
the resonance structures and the VER rate of the protein. 



4 Concluding Remarks 

In this paper, we have described theoretical (Sec.|2) approaches to the study of VER in proteins. 
We have examined VER of a CD stretching bond (mode) in cytochrome c from the QCF 
approach (Sec. 12.1^ and the reduced model approach fSec. I2.l3|) . For the CD mode in cyt 
c (in vacuum) at room temperature, both approaches yield similar results for the VER rate, 
which is also very similar to an estimate derived from an experiment by Romesberg's group. 
Our work demonstrates both the feasibility and accuracy of a number of theoretical approaches 
to estimate VER rates of selected modes in proteins. 

There are advantages and disadvantages of the (a) QCF approach and (b) reduced model 
approach to the prediction of VER rates in proteins. The QCF method is simple and applicable 
even for a large molecule like a protein. However, the VER mechanism may not be known a 
priori, and it must be supplemented by other methods such as anharmonic coefficient calcu- 
lations. Furthermore, the method relies on the local mode picture, which is easily applicable 
for high frequency (localized) modes, but not for low frequency (delocalized) modes. The re- 
duced model approach is quantum mechanically exact, and easily applicable for VER of low 
frequency modes. However, the anharmonic coefficient calculation is cumbersome even for the 
third order coupling terms in cyt c. Moreover, such a Talyor series expansion has not been 
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that the classical VER dynamics using an isolated methionine does not seem to be affected by 
including the fourth-order coupling elements (see the left of Fig. but we need to examine 
this issue further with quantum mechanical approaches such as TD-SCF methods. There is 
also an unsolved problem of the width parameter. Actually, this problem is not peculiar to 
the reduced model approach. The introduction of the width corresponds to coarse-graining, 
which also appears in the QCF approach when one averages the power spectrum of the force 
auto-correlation function. The most "ab initio" approach to solve this problem is a rigorous 
quantum mechanical treatment of the tier structure of energy levels in the protein |54j . The 
other appealing way is to regard 7 as a hopping rate between conformational substates j55| l5Hj . 
or a frequency correlation time, that may be derived from estimates of the frequency fluctuation 

[saisH]. 

Since both the QCF and reduced model approaches are based on Fermi's golden rule, there is 
a limitation for the strength of the interaction between the system and bath. There is a need to 
develop other methods without this deficiency. Promising approaches include nonequilibrium 
molecular dynamics methods |59j . time-dependent self-consistent field methods j32j . mixed 
quantum-classical methods |6U| . and semiclassical methods |27| . Another important issue is to 
calculate not only VER rates, but the physical observables related to the experiment data such 
as absorption spectrum or 2D-IR signals :61, 62 . In this case, we also need to deal with the 
effects of dephasing (decoherence) as well as VER. 

The accuracy of the force field parameters is the most annoying problem. Our preliminary 
calculations show that the VER rate in cyt c can vary by two order of magnitude when we 
change the bond force constant by ten percent (see the right of Fig. |3J. This situation is rather 
similar to that of the reaction rate calculation, where one must determine the activation energy 
accurately. Any inaccuracy in the activation energy causes an exponentially large deviation in 
the rate constant. This problem will be solved through ab initio quantum dynamics fSec. l2.2T^ 
or the reparametrization of the force field using experimental data or accurate ab initio calcula- 
tions. Given sufficient accuracy in the force field, we will be in a position to discuss the relation 
between the VER and function of a protein such as electron transfer in cyt c. As is well known, 
the dynamics of proteins related to function are well described by large amplitude (and low 
frequency) principal components [631 164j . The connection between principal components and 
VER should be investigated. The ergodic measure 65 will be a good device to examine this 
issue. As suggested by experiments |l4jn6j, collective motions in proteins can be important for 
the fast VER in proteins. The collective motions near the protein surface including solvation 
dynamics of water [661 1671 1691 l7Uj might be relevant for the VER and function. Cytochrome 
c and Mb remain excellent target proteins to investigate these fundamental issues of protein 
dynamics and its relation to function. 
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